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The Enskog-like kinetic approach, recently introduced by us to study strongly inhomogeneous flu- 
ids, is reconsidered in order to improve the description of the transport coefficients. The approach 
(T^ is based on a separation of the interaction between hydrodynamic and non-hydrodynamic parts. 

^-H The latter is treated within a simple relaxation approximation. We show that, by considering the 

C^ non-hydrodynamic part via a weighted density approximation, we obtain a better prediction of the 

CNj transport coefficients. By virtue of the simplicity of the kinetic equation we are able to solve numer- 

» I ically the phase space distribution in the presence of inhomogeneities, such as confining surfaces, 

^ I via a Lattice Boltzmann method. Analytical estimates of the importance of these corrections to the 

.^^ transport coefficients in bulk conditions is provided. Poiseuille flow of the hard-sphere fluid confined 

between two parallel smooth walls is studied and their pore-averaged properties are determined. 
O 

m 

I I I. INTRODUCTION 

^ The technological and industrial relevance of processes where molecular fluids flow through channels of nanometric 

«H diameter require to develop methods to control and predict such phenomena [T}{B]. Among the various techniques 

available today, one can mention continuum-based computational fluid dynamics methods, finite-element methods for 
Navier-Stokes equations ||7j , Molecular Dynamics [8, i9j , Kinetic theory LlQj Dynamic Boltzmann free-energy functional 
theory [TT], and Dynamical Density Functional theory for dense atomic liquids |12l [TB] . In these methods, one 
c/2 studies either the evolution of the local particle density n(r, t) or the evolution of the one-particle phase-space density 
distribution function. 

The strategy based on the solution of the phase-space distribution function is known as the Lattice Boltzmann 

method |14l I15| and is based on a discretized version of the Boltzmann transport equation, originally conceived to 

solve at a reduced computational cost the fluid dynamics equations. It can deal with complicated geometries and is 

^ numerically robust but, conversely, is typically based on top-down strategy: the phenomenological parameters, such 

Q as the viscosity and the thermal conductivity, are utilized as lumped parameter whose values are provided by some 

O experiment or macroscopic model. Since the underlying microscopic level is unexplored, in principle one cannot infer 

the nature of the atomic constituents. In addition, since the phenomenological transport coefficients are fixed from 

^^ bulk measurements, it is awkward to control how they change under inhomogeneous conditions, such as those due to 

^ confinement in narrow spaces. 

^^ In a series of papers we have considered an alternative approach [16-T8J, that does not require the a priori knowledge 

OO of the transport coefficients, nor the equation of state, but directly tackles the Boltzmann equation corresponding 

to a well defined microscopic model. A convenient route is to consider a system of hard-spheres and possibily 

adding attractive potential tails. The Hamiltonian determines the cross-section entering the Boltzmann-Enskog (BE) 

^T equation, thus establishing a link between microscopic and macroscopic properties. Using such bottom-up approach, it 

^^ is possible to determine all transport coefficients and the equation of state from the BE equation, without introducing 

, empirical ingredients. Such a strategy was pioneered by Luo, who carried out a systematic derivation of the Lattice 

. . Boltzmann equation starting from the Enskog equation. He obtained the equation of state for non-ideal gases together 

K** with thermodynamic consistency |19| . 

k>( The present method not only takes into account the non-local nature of the collision operator and leads to a 

"V^ non-ideal gas equation of state, but also reproduces the dependence of the transport coefficients on the density. 

C^ The approximation scheme utilizes concepts inspired by the Density Functional Theory in the weighted density 
approximation ( WDA) version |20l I21| , as regards to the use of a smoothed density field and the equilibrium pair 
distribution function as key quantities. At equilibrium our approach gives the same equations of the WDA and for 
such a reason we name it Weighted Density Lattice Boltzmann method (WDLB). Finally, the method is employed to 
derive a practical numerical scheme to study liquids fiowing in confined geometries. 

The purpose of this article is also to improve the WDLB by considering corrections to the viscosity and to the 
thermal conductivity. Such corrections are already present in the Enskog theory and, since they arise from the 
coupling between kinetic and coUisional modes, they are named cross-over contributions to the transport coefficients 
and depend linearly with the fiuid density. To this aim, we follow the clever approximation introduced by Dufty, Santos 
and Brey (DSB) ^22j to simplify the revised Enskog theory (RET) equation [23J, without spoiling its good features. 
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In particular, this approach preserves the non-locality of interactions, it accounts for the static pair correlations in 
a realistic way and provides a reasonable equation of state. In essence, Brey and coworkers proposed to replace the 
complex collision kernel featuring in the original RET (see eq. p4| below) by a more tractable object. 

The idea behind this proposal is to separate the relevant slow hydrodynamic modes from the fast evolving and less 
relevant kinetic modes. The first ones are then carefully handled, whereas the latter are treated by resorting to a 
simple relaxation time approximation, thus avoiding the effort of computing difficult integrals involving high velocity 
momenta of the collision operator. The resulting approximated collision operator consists of two pieces, one acting on 
the hydrodynamic modes and the second determining the relaxation towards equilibrium of the kinetic modes through 
the simplest approximation of the Boltzmann equation, the Bhatnagar-Gross-Krook (BGK) kernel |24| . 

Through this simplification, one obtains a representation where a set of effective fields describes the effect of 
intermolecular interactions in a self-consistent fashion, and the resulting transport equation is amenable to numerical 
solution. In a series of papers [25 27J we have numerically implemented this approach [22 by developing a numerical 
algorithm in the framework of the Lattice Boltzmann method. Here, we consider the effect of a more refined treatment 
of the relaxation term, as proposed by Santos, Montanero, Dufty and Brey (SMDB) [28J in order to account for the 
cross-over contributions to the transport coefficients. 

This paper is organized as follows: in section ITT] we briefiy review the employed kinetic method in order to present 
the extension of the WDLB approach, and motivate the introduction of the new terms into the equation of evolution 
used in numerical applications. In section [TTT] we give the explicit representation of the collision operator in the case of 
hard spheres and compute the effective fields. In section [TV| we present results for the transport coefficients. Finally, 
in section [V] we present the numerical test of the new approximation. In section [VT| we make some concluding remarks. 



II. MODEL 



The evolution of the phase-space one-particle distribution /(r, v, t) is represented by the following transport equation 



5t/(r,v,i) + v,dj{r,v,t) 



m dvj 



/(r,v,t)=17[/](r,v,t) 



(1) 



where F is an external force and $![/] represents the effect of the interactions among the fluid particles which we shall 
specify later. We adopt the Einstein summation convention that repeated indices are implicitly summed over and the 
notation di to indicate the partial derivative w.r.t. the i-th component of the vector r and dt to indicate the partial 
derivative w.r.t. time. 

In principle, after giving an explicit representation of 17 according to the underlying microscopic model, equation 
(fll can be solved numerically. The case of non-local functionals fJ, such as the RET, is computationally demanding. 
Alternatively, one can consider very simple forms of il, the BGK being the simplest ad hoc choice. Our strategy 
considers a non-local functional fi, but factorizes the velocity dependence in terms of a local Maxwellian distribution 
times Hermite tensor polynomials [25]. In addition, since one is chiefly interested in the evolution of the hydrodynamic 
variables it is possible to further simplify the form of Q, without altering the local transfer of momentum and energy. 

In this framework we implicitly assume that after few molecular collisions the system reaches a state of local 
thermodynamic equilibrium characterized by the hydrodynamic variables ri(r, t),u(r, i),T(r,i), being the number 
density, local fluid velocity and local temperature respectively. We further assume that the phase space distribution 
depends on space and time only through these flelds and their gradients. 

The flve hydrodynamic flelds can be obtained from /(r,v,t) through the relations 



n(r, i) 

n(r, i)u(r,i) 

|fcBn(r,t)T(r,i) 

where ks is the Boltzmann constant. 

We also introduce the corresponding moments of the collision operator 




/(r,v,t). 



(2) 



C(r.f) \ = jdw\ m- 

Q{r,t) J J \^m(v^ 



-u(r,t))" 
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r!(r,v,i). 



By multiplying both sides of eq. 
and integrating w.r.t. v we obtain t 



^ by the lowest Hermite tensorial polynomials (that is, 1, v, and m(v — u) 
le following set of balance equations for the flve hydrodynamic variables: 



(3) 



72) 



and 



dtn{Y, t) + 9,(n(r, t)u^{Y, i)) - , (4) 

mn{v,t)[dtUj{v,t) + Ui{v,t)diUj{v,t)\ + diPij{r,t) - Fj{r)n{r,t) -C-,(r,i) =0, (5) 

^n{r,t)kB[dtT{r,t) + u,{r,t)d.,T{r,t)] + P,,{r,t)d,uj{r,t) + dM^,t) - Q{r,t) = (6) 

The inoincntum and temperature equations dsl and ([6]) contain two additional quantities, a pressure term 

Pjj(r,i) =m / dv/(r,v,i)(v-u),(v-u)j (7) 

and a heat flux term: 

Ur, t) = 'jl '^v/lr, V, t)(v - u)2(v - u), (8) 

which in the n — >• Hmit can be identified with the kinetic components of the pressure tensor and heat flux, respectively. 
As we shall see in section |III[ these two quantities also contain contributions which are neither purely kinetic nor 
purely coUisional. 

One can also identify the last term in equation (IS]) with the gradient of the coUisional contribution to the pressure 
tensor, Pj^, 

air,t) = -l5,i^,(r,i) (9) 

and Q with a combination of the gradient of the coUisional contribution to the heat flux vector g^ and the pressure 
tensor times the strain rate: 



Q{r,t) = ~-[d,q'j{r,t) + P^,{r,t)dM^,t)]. 



(10) 



The second equality in eqs. ^ and (10 1 identifies the effective fields Ci and B with the thermodynamic fields. 

Formally, we have derived a set of equations for the five hydrodynamic variables, but these do not form a closed 
system, because we still need either to specify the non-equilibrium distribution function / or to assume some phe- 
nomenological constitutive relations between the gradients of Pij and q (the sums of the kinetic and coUisional parts) 
and the fluxes. Using the latter strategy one obtains the Navier-Stokes equations |32| . 

The alternative route is to focus on the evolution of /, with the collision operator fi being in general a non-linear 
functional of /, so that any kind of analytic work results hard. The simplest and crudest way to proceed is to 
approximate il by the following relaxation ansatz as done by BGK |211 

f7^«^(r,v,t)--^o[/(r,v,t)-n(r,t)0M(r,v,i)], (11) 

where 

, , , , m ,o/, / m(v — u)^ 

'^-(^'-'^)-[2^fc.7Xr;7)]'^^H"2Wr;^ 

is the local Maxwellian. The BGK recipe for fl preserves the number of particles, the momentum and the kinetic 
energy, in other words it fulflUs the physical symmetries and conservation laws of the fluid. Finally it gives a quite 
uninteresting thermodynamic behavior corresponding to an ideal gas. Instead, it gives non vanishing values of the 
transport coefficients, which turn out to be functions of phenomenological parameter wq |24| . 

In view of the interest for the non-ideal gas behavior one seek for a better approximation capable of predicting 
a non-trivial equation of state and transport coefficients. The Boltzmann equation for hard spheres has non-trivial 
transport coefficients but lacks a satisfactory equation of state. The RET form of O has the sought features, but is 
hardly tractable in numerical work [23j. 



An interesting and practical reduction of the RET collision operator has been devised by Dufty et al. [22] . Their 
idea is to separate the contributions of fi to the hydrodynamic equations from those affecting the evolution of the 
non-hydrodynaniic modes, by projecting the collision term onto the hydrodynamic subspace spanned by the functions 
{l,v, u^} and onto the complementary kinetic subspace: 



with 



^' ' hydro^^ i V-* ' hydro)^' 



-Pkydron = i^^^^^M (r, t) [(V - U) • C(r, t) + [ "^(^^^"(^^^))' _ l]Q(r, t) 



The second term, il± = (/ — Vhy 



SkeT 
:,)fl is further approximated as: 

n± sa -a;o[/(r, v, t) - n(r, t)(j)M{r, v, t)] 



(12) 



(13) 



(14) 



and (ITl by a simpler equation, where the complicated interaction between non hydrodynamic modes is approximated 
via a BGK-like relaxation term with ujq being a phenomenological collision frequency chosen as to reproduce the 
kinetic contribution to the viscosity. 

It was pointed out by Santos et al. |28^ that better results for the transport properties could be obtained upon 
modifying the BGK relaxation term in ( 14 1 as to allow the relaxation time to depend on the spatial inhomogeneities 



of the local equilibrium state. This effect depends on the coUisional transfer mechanism associated with the difference 
in position of the colliding particles. To achieve this goal, following Santos, Montanero, Dufty and Brey [2S], we 
complement the relaxation term in eq. (14 1 in the following way: 



f7i*^^^(r,i) = -c^o[/(r,v,i)-n(r,i)</'M(r,v,i)] + (m^)^0M(r,v,i)x 



(v-u) 
r?T,/3(v - 



? 5 



3(v-u)25,, 



;Ay(r,t) + 



(v-u),^S,(r,t)} 



(15) 



The extra terms in eq. ( 15 1 give a better representation of the decays towards equilibrium as discussed in ref. 



The first term is proportional to ujq, which assumes a constant value throughout the system, whereas the second term 
may depend on the spatial inhomogeneities of the velocity and of the temperature fields. Such a dependence occurs 
through the functions Aij{r,t) and Bi{r,t) and thus is modulated by the spatial structure of the fluid. In the case 
of a uniform system the relaxation process may occur only via the term proportional to wo, when the distribution 
function / is different from the local equilibrium distribution. The inclusion of these terms also yields a quantitative 
improvement of the transport coefficients with respect the simpler theory discussed previously. 
The new terms can be obtained by computing the following projections of the collision operator: 



and 



^„-(r,<) = / dv[(v - u),(v - u)j - glv - ufS,j] f7(v,r,i) 

B,ir,t)^ /'dv[(v-u)2-5^](v-u),;f^(v,r,i) 
J m 



(16) 



(17) 



By construction, the new terms do not affect the structure of the continuity equation m, momentum balance ([5| 
and energy balance equations (l6|. The SMDB equation can be cast in the form: 



dtf{r,^r,t) + v,dif{r,-v,t) 



F^ir) d 



f{r,v,t) = nr^-ir,t) 



+ m/?</>A./(r,v,t){(v-u)-C(r,t) + ( 



m/3(y 



1 



)Q{r,t)} 



(18) 



where /3 = l/kBT{r, t). A simple calculation allows to compute the deviation 5f of the distribution function /(r, v, t) 
form its local equilibrium value fioc = n{r,t)(f>M{'i^,t). We replace / — >• fioc in the left hand side of eq. (18 1 and 



obtain: 



I \dtn + di{nui}] + \ndtUk + nuidi{nuk) + dk{nT) Ck 



Fk „ 1 rn{vk - Uk) 



knT 



rn/? / (v — u)'^ \ ~l 
+mP[ndiUk —Aik]\{v., ~ Ui){vk - Ufc) %j j^Af (r, v,i) = -WQi^/Cr, v, t) (19) 

Since the projections of 6f over (1, {vi — Ui), m(v — u)^/(A;bT) — 3)) vanish , we impose that also the first three terms 
in the l.h.s. vanish. These are the so-called solvability conditions and are precisely the balance equations (Bl, (Is]) and 
(l6| . We thus obtain the following explicit representation of 5f 



-I- m/3 n{v,t)diUk{r,t) T^Aik \ {vi - Ui){vk - Uk) 5ik]\ (20) 



2 kBT{v,t) 

)} 



The fields Aij and Bi represent a measure of the distortion of /(r, v, t) with respect to the local equilibrium fiod^^, v, t). 
Using such an approximation for df, we can now estimate how the distortion affects the components of the pressure 
tensor (|7|: 

P,,(r,i) = kBT{r,t)n{r,t)[d,j - ^(^{d,Uj{r,t) + djU,{r,t)) - ^% ^^^^^(r, t)) }+(5Py(r, t) (21) 



where the first term is the BGK result and 

m 

OJQ 



cy 

5P^j{v,t) = — A,j{Y,t)+Aj,{Y,t) - -%^Afcfe(r,t) (22) 
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stems from the Aij contribution. The heat flux is obtained by substituting 5f in eq.(|8| 

ft(r, t) = -l-^n{Y, t)k%T{r, t)d,T{v, t) + 5q,{r, t) (23) 

with<5g,(r,t) = ^fi?,. 

It is worth noting that Luo [33] compared the approaches derived from the Enskog equation with those based on the 
simplifying assumption that the intermolecular forces are assimilated to a forcing self-consistent term. He concluded 
that the latter are plagued by an inconsistency. In fact, the pressure featuring in the momentum equation must also 
feature in the energy equation, a condition which is not satisfied if one treats the interactions in a mean-field fashion 
[Mll35| . In this respect the present treatment has the correct form, since eq. (l6| contains the viscous dissipation term 
representing the irreversible transformation of kinetic energy into heat. 

III. HARD SPHERE MODEL 

At this stage we are interested in obtaining an explicit representation of the fields {Ci,Q,Aij, Bi) in order to 
determine the viscosity and thermal conductivity, and to implement the numerical solution of the evolution equation 



(18 1 by lattice methods. To this purpose we use as a benchmark the prototypical model of dense fluids, namely the 
hard-sphere system at moderate packing fractions. The RET collision operator involves only the single-particle phase 
space distribution and accounts for static pair correlations, but not velocity correlations. At high density it becomes 
more inaccurate because it disregards the occurrence of correlated collisions. Nevertheless, the RET approach is very 
accurate for our purposes and has the interesting feature of taking into account the instantaneous transfer of energy 
and momentum between particles whose distance equals the hard sphere diameter. For this reason even to a spatially 
uniform state corresponds a non-ideal gas pressure. 



The RET collision operator reads j23) : 

r!^^^(r,vi,f) = a'^-i /'dv2 /"dfce(k-vi2)(k-vi2)x 

{52(r,i 



,r-ka,t)/(r,v'i,t)/(r-ka,v^,t)- 
52(r,r + kfT,i)/(r,vi,i)/(r + ka,V2,i)} (24) 



where v[ and Vj are scattered velocities determined from v'l = Vi — (k • Vi2)k and Vj — V2 + (k • Vi2)k, a is the 
hard-sphere diameter, k is the unit vector directed from one particle to another. Q is the Heaviside function, and 
.92 (ri, r2, i) is the pair correlation function evaluated at contact distance |r — r'| = a. It describes the positional pair 
correlations corresponding to the inhoniogeneous single particle density profile n{r, t). There is no velocity dependence 
in the two particle correlation. 

Following Ref. [33] we assume that (72 at time t is the same as the one corresponding to the same system having 
a density profile equal to n(r, t) at thermodynamic equilibrium, an assumption often referred to as the adiabatic 
approximation. Moreover, since the exact inhoniogeneous form of 172 is unknown, we resort to the Fischer-Methfessel 
prescription [37], stating that g2(ri,r2ji) is given by the pair correlation function at contact of an homogeneous 
system evaluated replacing n by the following smeared density: 



n{r,t) = — - / dr'n(r + r'). 

•^O" J\r'\<cr/2 



(25) 



Because we have access only to the first lower velocity moments of /, we cannot solve the full problem (IT]) using 
such a collision operator. The DSB and SMDB recipes represent viable options. Recently we have used the SDB 
equation in conjunction with the Lattice Boltzmann method to study one component, two component and ternary 
mixtures with interesting applications also to electrolyte solutions [38j ,39^. Hereafter, we want to test the SMDB 
equation implemented by an appropriate LB method and compare the results. 

The expressions of C, i? have been obtained before [TT and read 



and 



(26) 



Ci{r,t) — ~-kBT{r,t)n{r,t)(7 / dkkig2{r,r + a'k,t\n)n{r + akjt) 



^7rfcBT(r,i)/m 



k ■ [u(r + crk,t) -u(r,i)] + 



T{r + ak,t)-~T{r,t) 
2r(r,t) 



Q{r,t) ^ kBTir,t)n{r,t)a'^ / dkg2{r,r + a'k,t\n)n{r + a'k,t) 



fc- [u(r + fc,t) -u(r,i)] 1 kBT{r,t) T{r + ak,t) - T{r,t) 
2 ^T^V m T(r,i) 



(27) 



A similar calculation was performed to evaluate Aij and Bi using (16 1 and (17|, respectively, under the approxi- 
mation that n^^'^lf, f] featuring in these two equations was replaced by fi-'^-^-^f 



Xcfioc] mM- We find: 



A,,(r,t) 



kBT{r, t) 
5m 



n{r,t)a^ / dkg2ir,r + ah,t\n)n{r + ah,t) 



kj[ui{r + crk, t) — Ui{r,t)] + ki[uj{r + ak^t) — Uj{r,t)] 



Similarly, we derive the correction to the heat fiux vector 



B, (r,t) = -2-%r(r,t)7i(r, i)cr^ / dkkig2ir,r + ak,t\n)n{r + ak,t) 



T{r + ak,t)^T{r,t) 



(28) 



(29) 



As shown in ref. [26 a useful representation of the term (26 1 is obtained by decomposing it in the following way: 
C{r,t) — n{r,t) (F"^^{r,t) + F"*'"^°"*(r,i)), where we can identify a force acting on the particle at position r as due 
to the influence of all remaining particles in the system, the so called potential of mean force ^41] : 



F"^(r,t) = -fesTcr^ / dkkig2{r,r + ak,t)n{r + ak,t) 



and a viscous term 



2 rnksT 



dkkikjg2{v. r + uk, t)n(r + crk, t)(uj(r + crk) — uj{y)) 



(30) 



(31) 



and a force due to the presence of thermal gradients 



Ff{Y,t) = -^ I dkk,g2{r,r + ak,t)n{r + ak,t)kB[Tir + ak,t)-T{r,t)] 



(32) 



One can see that the form of the fields Ci,Q,Aij,Bi depends on density, velocity profile and temperature. Only 
one contribution to Ci inherits the property of the density functional theory of being fully determined by the density 
only. The other terms require the knowledge of the hydrodynamic fields even in the present approximation, since they 
are not slaved to the density field as it occurs in Dynamical Density Functional theory. 

IV. ANALYTIC ESTIMATE OF THE BULK TRANSPORT COEFFICIENTS FROM WDLB 

Let us consider the limit of small velocity gradients and constant density and temperature, so that: 



A,{r,t)^-'^n^g,ia+)^a' 



and, using eq.(22l: 



m 



1 • -.2„ . +N87r , 



. 47r 
— ( 
15 



djUi + diUj 



5P^,{v,t) = kBTn'g2i<J+)-^<7 



ujo 



15 



2 .r-^ 

djUi + diUj - -5ij 2_^ dkUk 



(33) 



(34) 



Let us consider a system of uniform density and temperature, and a small sinusoidal velocity perturbation varying 
in a direction normal to the streamlines: 



By neglecting the non-linear term, we obtain from eq. ([5| 

d d ~ 

— Uy{x,t) = -—Pxy{x,t) +Cy{x,t) 

and, substituting the expressions for Cy and P^y, 



n—uy{x,t) 



r 1 



Stt 



kBTn \ 1 + -—na^g2{a^) ) + ■— \J rm:kBTna^ g2{<7^) 
Wg \ 15 / 15 



d'^Uy{x,t) 
dx'^ 



(35) 



(36) 



(37) 



The solutions have exponential form, Uy{x,t) — Uj,(0,i)e'^^^e^'^''/"^''=*, where the shear viscosity 77 is given by the 
expression inside square brackets in the r.h.s. of eq. pTJ). 

The first term arises from the reduction of the velocity gradients due to the displacement of particles from a region 
characterized by a higher fluid velocity u to a region of lower velocity and vice versa. The third term takes into 
account the fact that neighboring volume elements moving at different velocities can exchange momentum due to the 
interaction of their surface molecules without involving particle displacement. The second term 77^^), cal led cross- over 
term, has a mixed nature |32]- Finally, we set the collision frequency to the value ujq = ^nu^g2 (ct) yJi^kBT jra |28| 
giving the correct low density shear viscosity of hard spheres of diameter a 



^ ^ r^^ ^^ ^ 8j^^3^^^^+^^_^ 4 ^^^^^^^--^^^4^^^^^ 



Wo 



15 



15 



(38) 



A. Bulk viscosity 



In order to derive the expression for the bulk viscosity, 77^, we consider a velocity field of the form u — {ux{x, t), 0, 0), 
use the macroscopic relation 



dP^Ax,t) , , 4 ,a2u,(x,i) 



to define 77^ and obtain after a simple calculation 



dPxxix,t) 
dx 



4 1 



- — kBTn [ 1 + -—na^g2{a^) ) + -^/TT\/rnk^na^g2{a^) 
3 Wo \ 15 / 5 



d'^Ux{x,t) 
dx'^ 



(39) 



(40) 



By comparing, the two expressions for the derivative of the pressure tensor, we find 7/t, ~ ^^77 and conclude that rjf, is 
purely coUisional with no kinetic and cross contributions. 



B. Thermal conductivity 

The heat flux correction Bi for small gradients can be approximated as 

SUr,t) c -^'^T{r,tWg,{a+)^a^^ (41) 

ujQ m b ox 

Let us consider the heat flux in a system at rest and of uniform density. According to equation m 



fcBn(r,t)atr(r,0 = {^^n(r,t)4T[l + ng2(a+)^a3]+^y;^^;^fc^52(a+)nV^}|J^ (42) 



yielding the thermal conductivity, A, as the sum of the three terms inside the parenthesis in the r.h.s. of (42 1. These 
are the kinetic, cross-over and coUisional contributions to A, respectively. For the sake of comparison the Enskog 
formulae [42_ are 

riE = 77^") [1.016/32 (cr+) + 0-Sbn + 0.761g2{a^)b'^n^] (43) 

and 

A_E = A(°' [1.025/32 (cr+) + 1.230 bn + 0.776Ag2{a+)b^n'^] (44) 



where 6 = 27r/3a3 and 77(0) = ^J^^^^i;^ and A(o) = ^^ ' ''"^ 



We shall derive such corrections and compute numerically their effect by using the SMDB approach. 



NUMERICAL RESULTS 



By using the results derived above, with the effective fields given by (26), (27 1, (28 1 and (29 1, we construct the 
corresponding Lattice Boltzmann algorithm to solve the kinetic equation ( 18 ) on a three-dimensional cubic lattice 
with lattice constant a, assumed to be a fraction of the hard sphere diameter a. 

The continuous dependence of /(r, v,i) on its argument r is replaced by a finite set of nodes of coordinates r^ of 
a regular lattice spanning the physical space available to the particles. The v dependence is given by spanning / 
over an orthonormal basis formed by a finite number of tensorial Hermite polynomials. We consider the so-called 
D3Q19 version, where the v-space is characterized by 19 discrete velocities. The state of the system is specified by 
the discrete set of distribution functions fa{ri,t) and for the sake of simplicity, we only consider isothermal systems 
and neglect the evolution of the T field. The finer details of the method have been described in refs. |16} ITT] and will 
not be repeated here. In the present implementation, the method only contains the two additional terms proportional 
to Aij and Bi. 

The LB algorithm performs the following basic operations for every mesh node: 

1. the discretized distribution is initialized by specifying the local values of the density and fluid velocity; 



2. the values of the flelds Aij ,Bi,Ci,Q are computed via eqs. ( 26 )-( 29 1 using the fields u and n; 

3. the distribution is evolved in time though the collision and streaming steps via an explicit trapezoidal rule; 

4. the new fields n(r,t) and u(r,i) are computed using the evolved distribution function; 

5. the cycle is repeated starting from step (2). 

In this scheme, one relevant question pertains the choice of the BGK relaxation frequency wq. According to the Enskog 
picture, the kinetic contribution to viscosity should be set equal to 7]k — ji^pixl ^^'^ and therefore wp = nhsX^ j^ 
addition, the low-density regime should be characterized by a constant dynamic viscosity, as first observed by Maxwell 
and, as n — 7> 0, the viscosity should go to zero. This rather large baseline for viscosity necessarily implies a strong 
departure from the local equilibrium. Such a condition is numerically critical, as it corresponds to a mean free path of 
the order of a, whereas numerical stability imposes that the mean free path should be of the order of the mesh spacing 
a. In this case, in fact, the BGK term becomes effective in tethering the distribution close to the local equilibrium 
form, coherently with the small-Knudsen conditions employed in the theory developed so far. For this reason, we 
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choose a smaller value of the kinetic viscosity. It is important to remember that also the cross-over term depends on 
wq, but nevertheless we can still determine its effect on the global viscosity. 

In the numerical benchmark, we evaluate the effect of the cross-over contributions to shear viscosity stemming from 

I 3 

eq. (15 1. In this test we choose a hard-sphere diameter a = 8a and the packing fraction x — ^'^g " ranges from dilute 
gas to dense liquid conditions (% ~ 0.40). In Fig. fTlwe report the shear viscosity in bulk conditions as a function of 
the packing fraction x- The measurements were performed by initially perturbing a uniform quiescient state by a sine 



wave of the type u^ sin^q^x) and monitored that its relaxation occurred with a characteristic time T{qx) = ~Qxj where 
r] is the shear viscosity and hq the density of the uniform fluid. The data are overall smooth up the highest packing 
fraction and we observe that the correction to viscosity is relevant. The presence of the crossover term determines a 
larger viscosity over the entire range of packing fractions studied with respect to the one obtained in the absence of 
it. 
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FIG. 1: Shear viscosity computed in bulk conditions ribuik as a function of the packing fraction (this quantity should not be 
confused with the bulk viscosity of Section 4.1). Circles and squares are with and without the cross-over correction to viscosity, 
respectively. 



Next, we address the classical question regarding how the confinement in a narrow space affects the viscosity of 
the fluid |151 - li5] . In our numerical tests, we consider the flow of a hard-sphere fluid between two parallel plates 
and observe the density and velocity profile. Since this a direct test of the viscosity of the fiuid we expect a change 
with respect to the previous version where the terms Aij and Bi were not taken into account. As shown in FigJ2] 
the measured viscosity converges for quite large channels of the order of 200 lattice units towards the predicted bulk 
value, whereas for small values of the width W is consistently lower, thus indicating that the walls strongly affect the 
motion of the adjacent layers. 



VI. CONCLUSIONS 



In this study we have considered the corrections to the transport coefficients which determine the cross terms 
between kinetic and collisional contributions. At the price of including a couple of extra contributions in the term 
describing the relaxation of fast modes, that is by lifting the BGK approximation for these modes, we have obtained 
a new estimate for the transport properties. After deriving an explicit expression for these quantities, we have 
evaluated the correction to viscosity via numerical simulations. It is interesting to remark that the range of validity of 
the presented LB algorithm does not extend to very low densities. The reason being that, if one assumes a constant 
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FIG. 2: Normalized shear viscosity as a function of the effective slab width {W 
squares correspond to packing fractions of 0.11 and 0.32, respectively. 



a) expressed in lattice units. Circles and 



relaxation frequency wq , the value of the kinematic viscosity is too small and the small Knudsen number assumption 
breaks down. 

The resulting form of the expression for the viscosity has similarities with the one derived within the local average 
density model (LADM) by Davis and coworkers [46 . They assumed that the functional dependence of transport 
coefficients on density can be computed via the corresponding expressions in the Enskog theory of hard-sphere fluids, 
by replacing the actual density n(r,t) by the coarse grained density, n{r,t). An approximation very similar to the 
LADM can be derived within our method by using the explicit form of Q. Such a recipe is in the same spirit as 
the early weighted density density functional theories of equilibrium systems where the local free energy density was 
assumed to be a function of n(r). Moreover, as pointed out by Hoang and Galliero in a Molecular Dynamics study 
|47j . the LADM underestimates the importance of the local density as far as the kinetic term is concerned. The 
present theory, does not suffer from such a problem since it treats on a separate basis kinetic and potential terms. 

Finally, the presence of attractive interactions can be accommodated in a mean field fashion by adding a Random 
Phase term in the evolution equation. However, while representing a very important contribution to the equation of 
state, such a modification does not modify the transport coefficients (with the notable exception of the interspecies 
diffusion coefficient in mixtures, see refs. |27ll48j V The reason for this shortcoming is that in the Enskog-like equation 
the details of the potential enter only through the cross-section. 
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